

function tfr_ISI(sb)

gavgdir = '/home/predatt/anatod/ISI/MEG/GAVG';
cd(gavgdir)
load subjname
load dataset

rootdir = '/home/predatt/anatod/ISI';
datadir = '/home/predatt/anatod/ISI/raw'; 

infofile = strcat(subjname{sb},'_datainfo.mat'); 
subdir = fullfile(rootdir,'MEG',subjname{sb}); % subject directory
cd(subdir);
load(infofile)
cd(datadir)
    

% select triggers and conditions of interest	
cond_trig{1} = [33]; % focal, attended 1st
cond_trig{2} = [34]; % focal, attended 2nd
cond_trig{3} = [35]; % distribued, attended 1st
cond_trig{4} = [36]; % distributed, attended 2nd

for i=1:length(cond_trig)
        cfg = [];
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl;
        trial_indx = [];
        for j=1:size(cfg.trl,1)
            if sum(cfg.trl(j,4)==cond_trig{i})
                trial_indx(end+1)=j;
            end
        end
        
        index{i} = trial_indx;
   end
   
selection = [index{1} index{2} index{3} index{4}];

% cut out larger trials for better low frequency estimation        
        cfg = [];
        datainfo.trl(:,1) = datainfo.trl(:,1) - 1800; % first sample is shifted 1.5 sec
        datainfo.trl(:,2) = datainfo.trl(:,2) + 1800; % last sample is shifted 1.5 sec
        datainfo.trl(:,3) = datainfo.trl(:,3) - 1800; % offset to stim onset is shifted  1.5s
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl(selection,:);
    
   % preprocess 
        cfg = [];   
        cfg = datainfo.basicinfo;
        cfg.trl = datainfo.trl;
        cfg.continuous = 'yes';
        cfg.channel = {'MEG'};
        tmpdata = ft_preprocessing(cfg);        
        
        if sb == 17
            grad = tmpdata.grad;
        end
        
   % find neighbours
        cfg = [];
        cfg.method = 'distance';
        neighbours = ft_prepare_neighbours(cfg,tmpdata);
    
    % quickly redo ICA
        cfg = [];
        cfg.unmixing     = datainfo.comp.topo;
        cfg.topolabel    = datainfo.comp.topolabel;
        comp_data        = ft_componentanalysis(cfg, tmpdata); 
  
    % remove unwanted components
        cfg = [];
        cfg.component = datainfo.comp.artifact; 
        data = ft_rejectcomponent(cfg, comp_data);
        clear comp_data
        
        % preprocess 
        cfg = []; 
        cfg.channel = {'MEG'};
        cfg.demean = 'yes';                      
        cfg.blcwindow = [-inf +inf];          
        cfg.detrend ='yes';                 
        cfg.dftfilter = 'yes';              % notch filter to filter out 50Hz
        cfg.lpfilter = 'no';                % low-pass filter
        cfg.continuous = 'yes';
        data = ft_preprocessing(cfg,data);
        
                
   if sb == 17
        
        cfg = [];
        cfg.method = 'template';
        cfg.layout = 'CTF275.lay';
        cfg.neighbours = ft_prepare_neighbours(cfg);
        cfg.missingchannel = {'MRF32'};
        cfg.method = 'nearest';
        [data] = ft_channelrepair(cfg, data);
        
   data.grad = grad;
   end

    % Calculate planar gradient of each trial
        cfg = [];
        cfg.planar = 'yes';
        cfg.planarmethod = 'sincos';
        cfg.trials = 'all';
        cfg.neighbours = neighbours;
        cfg.method = 'template';
        data_planar = ft_megplanar(cfg,data);
        clear data;

%% loop over conditions
% 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Do TFR calculations using Hanning taper for low frequencies
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        cfglofreq = [];
        cfglofreq.output      =   'pow';
        cfglofreq.channel     =   'MEG';
        cfglofreq.method      =   'mtmconvol';
        cfglofreq.keeptrials  =   'yes';
        cfglofreq.keeptapers  =   'no';
        cfglofreq.taper       =   'hanning';
     
 %         adaptive taper length;
         cfglofreq.foi         =   [5:1:35];
         cfglofreq.t_ftimwin   =   [3./ cfglofreq.foi]; % 3 cycles
         cfglofreq.toi         =   [-0.5:0.025:1];   % time of interest
         adaptive taper length;
        
        data_tfr            = ft_freqanalysis(cfglofreq,data_planar);
        tfr_lo              = ft_combineplanar([],data_tfr); 

        clear data_tfr;
        cd(subdir)
        TFR_lo_filename = strcat(subjname{sb},'_TFR_lo');
        save(TFR_lo_filename,'tfr_lo', '-v7.3');
             
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % Do TFR calculations using multitapers for high frequencies
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % fixed time and frequency smoothing
%         cfghifreq = [];
%         cfghifreq.output      =   'pow';
%         cfghifreq.channel     =   'MEG';
%         cfghifreq.method      =   'mtmconvol';
%         cfghifreq.pad         =   'maxperlen';
%         cfghifreq.keeptrials  =   'yes';
%         cfghifreq.keeptapers  =   'no';
%         
%         cfghifreq.foi         =   [30:3:120]; % frequencies must fit in the time window; hence 3 Hz steps
%         cfghifreq.t_ftimwin   =   0.2 * ones(1,length(cfghifreq.foi)); % 200 ms
%         cfghifreq.tapsmofrq   =   10 * ones(1,length(cfghifreq.foi)); % 10 Hz
%         cfghifreq.toi   =   [-0.5:0.025:1];   % time of interest
% 
%         data_tfr        =  ft_freqanalysis(cfghifreq,data_planar);
%         tfr_hi          =  ft_combineplanar([],data_tfr);   
%         clear data_tfr;
% 
%         cd(subdir)
%         
%         TFR_hi_filename = strcat(subjname{sb},'_TFR_hi');  
%         save(TFR_hi_filename,'tfr_hi', '-v7.3');


end



